function [Ss, rs] = find_spillovers_mto(Ss, Hs, rs, myDist, tol, max_It, trans)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha...
        xi xi1 ub lb pi_e H w theta Int_pi_w sigmae elasts Pr log_e_mu log_e_sigma phi sigma gamma rts r_bar fp...
        r_base base_rs w_mult shares w_bar tmp_shift
   
    % Fixed Point Iteration to find spillover

    pi_w = sum(myDist, 2:3);
    w_bar = dot(w, pi_w);  
   
    
    for it = 1:max_It
        % Spectral Algorithm to find market clearing prices

        x_0 = rs(1:end);
      
        x_1 = x_0 + (size(rs,1):-1:1)'/100; 
    
        f_0 = find_rents_mto(x_0, Ss,  Hs, myDist, trans);
        f_1 = find_rents_mto(x_1, Ss,  Hs, myDist, trans);

        max_iter = 100;
        for i = 1:max_iter 
            s = x_1(1:end) - x_0(1:end);
            y = f_1(1:end) - f_0(1:end);

            x_1 = x_0;
            f_1 = f_0;

            lambda = dot(s,y)/dot(y,y);

            step_size = 1;
            for iter_sub = 1:20

                x_0_t = x_1 - lambda * f_1 * step_size;
              
                if sum(x_0_t(1:end-1) > 0 ) == (length(x_0_t)-1)
                    break
                else
                    step_size = step_size * 0.9;
                end
            end
            x_0 = x_0_t;
         
            f_0 = find_rents_mto(x_0, Ss,  Hs, myDist, trans);

            diff = sqrt(dot(x_0 - x_1, x_0 - x_1));
            diff_f = sum(abs(f_0));

            if (diff < 1e-6) & (diff_f < 1e-3)
                rs = x_0;
                break
            elseif i == max_iter
                
                rs = [w_bar*0.7, w_bar*0.5, w_bar*0.2]';
                x_0 = rs(1:end);
                minf = @(x) sum((find_rents_mto(x, Ss,  Hs, myDist, trans)).^2);
                options = optimset('TolFun',1e-6,'TolX',1e-3);
                [minx, fval, exitflag] = fminsearch(minf, x_0(1:end), options );
                
                rs = minx; 
                if (fval > 1e-3 | isnan(fval)) & (it > 15 )
                    fval, minx
                    error("Did not find root")
                end
            end                
        end
        [~, tax, mean_wage] = find_rents_mto(rs, Ss,  Hs, myDist, trans);
        % Calculate Choice Probabilities
        [P_S, ~, wprimeS, ~] = choiceProb_ct_par(Ss, rs, tax, trans, mean_wage);
        P_S = sum(bsxfun(@times, P_S, reshape(pi_theta, [1,1,1,1,2])),5);
        popMass = bsxfun(@times, P_S, myDist);

        wMass = squeeze(sum(bsxfun(@times, wprimeS, popMass), 1:3));
        popMass = squeeze(sum(sum(popMass, 1:3),5));


        Ssi = wMass ./ popMass;
        Ssi = Ssi(:); 
        Ssi = Ssi.^xi1;
        d_S = abs(Ssi - Ss);

        % Test for convergence
        if sum(d_S) < 2 * tol
            break
        elseif it == max_It
            display(sum(d_S(:)))
            error("Did not converge in Spillover") 
        end
        Ss = Ssi;

    end
    
end